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•O ' Abstract 

^: 

f~^ ^ Two fundamental quantities for characterizing nonlinear wave phenomena in plasmas are 

V^ ' the spectral energy transfer associated with the energy redistribution between Fourier modes, 

("^ ' and the linear growth rate. It is shown how these quantities can be estimated simultaneously 

0^ , from dual-spacecraft data using Volterra series models. We consider magnetic field data gath- 

CTN ■ ered upstream the Earth's quasiparallel bow shock, in which Short Large Amplitude Magnetic 

C/3 ' Structures (SLAMS) supposedly play a leading role. The analysis attests the dynamic evolu- 

O , tion of the SLAMS and reveals an energy cascade toward high-frequency waves. These results 

C/3 ' put constraints on possible mechanisms for the shock front formation. 

>>: 
^: 

P-t- 1 Introduction 

>'■ 

■ ^ ■ In this paper, a Volterra series representation is used to describe the nonlinear evolution in time 

r^J I and in space of a fluctuating wave fleld. The basis for this approach is that plasmas can often be 

S ■ viewed as a causal nonlinear system (a "black box") that reacts to a given excitation by giving a 

response. By modeling the nonlinear transfer function associated with this system, deeper insight 

can be gained into the underlying physics. 

The analysis of the nonlinear transfer function is detailed here for the particular case where two- 
point measurements are available. First, we show how to model the dynamical response. Then, the 



physical interpretation of the model coefficients is given. Particular attention is paid to the linear 
growth rate, which expresses the linear instability of the wave field, and to the spectral energy 
transfer, which describes how the instabilities saturate through nonlinear wave interactions. 

We apply this method to magnetic field data gathered by the dual AMPTE satellites near 
the Earth's quasiparallel bow shock. This data set corresponds to a regime of quasi-stationary 
turbulence in a coUisionlcss plasma; it has received much interest in relation with the existence of 
Short Large Amphtude Magnetic Structures (SLAMS) [Schwartz et al., 1992]. It is shown how a 
transfer function analysis reveals the role played by these nonlinear structures. 

This paper is divided in three parts. The experimental context is described in section ^ 
Sections ^ to are devoted to data analysis aspects with a description of the model, the choice 
of its parameters, and its validation. Finally, in sections g and y, experimental data are analyzed 
and interpreted. 

2 The Experimental Context 

The magnetic field data of interest were gathered by the dual Active Magnetospheric Particle 
Tracer Explorers spacecraft (United Kingdom Satellite (AMPTE-UKS) and Ion Release Module 
(AMPTE-IRM)) on day 304 of 1984 just upstream the Earth's quasi- parallel bow shock. Several 
studies have already been devoted to this particular event [Schwartz and Burgess, 1991, Schwartz 
et al., 1992; Mann et al., 1994; Dudok de Wit and KrasnoseVskikh, 1995; Dudok de Wit et al., 
1995], which provides a paradigm for nonlinear effects in turbulence. The spacecraft were closely 
following each other on the same outbound orbit (with a separation of 5x — 144 km), depicted in 
Figure ^ 
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Figure 1: Configuration of the Earth's bow shock showing the magnetic field lines, the orientation 
of the solar wind, and the orbit of the spacecraft for the event considered in this paper. 



A distinctive feature of the studied region is the occurrence of SLAMS, which supposedly play 
a leading role in the shock front formation [Schwartz et al., 1992]. The shock wave is caused by the 
sudden deceleration of the supersonic solar wind at the encounter of the Earth's magnetosphere. 



The SLAMS grow out of low-frequency waves that propagate away from the shock front but are 
convected back toward the Earth by the solar wind [Thomsen et ai, 1990]. This steepening process 
is likely to result from an interaction with ion beams coming from the shock front [Scholer, 1993]. 
There are several open questions regarding the role played by SLAMS. Quasi-parallel shocks 
are currently viewed either as an entity [Winske et ai, 1990] or as a patchy transition zone made 
by a merging of SLAMS [Schwartz and Burgess, 1991]. The relationship between the SLAMS 
and the whistler wave packets that frequently occur at their leading edge is not well understood 
either, although there is numerical [Omidi and Winske, 1990] and experimental [Dudok de Wit and 
Krasnosel'skikh, 19951 evidence for a causal link between the two. 
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Figure 2: Excerpt of the time evolution of the magnetic field amplitude, as measured by the two 
spacecraft. A typical SLAMS appears at t —71 s. The precursor whistler wave appears at its 
trailing edge since the wave field is convected backward by the strong solar wind. 



An excerpt of the magnetometer data is shown in Figure 0. The trajectory of the spacecraft, 
the prevailing magnetic field and the average solar wind velocity {vsw = 370 km/s) are all parallel 
within a few degrees. This is an important point since it means that both spacecraft see the same 
structures, separated by a time interval of about 1 s. A comparative analysis should therefore 
reveal how the wave field, and in particular the SLAMS, evolve as they move from one spacecraft 
to the other. We do this by building a Volterra model that tries to predict the wave field of 
AMPTE-IRM using the data of AMPTE-UKS as input. 

Each spacecraft provides a data set which consists of the three components of the magnetic 
field, measured immediately upstream the shock front. For each component the number of samples 
is 4521; the data were sampled at a constant rate of 8 Hz after being low-pass filtered at 4 Hz. We 
have chosen to consider the three components as different ensembles, thereby artificially increasing 
the sample size by a factor of 3. The anisotropy of the wave field a priori does not justify such 
an approximation, but no significant differences were found between the model coefficients as 
estimated separately from each component. An obvious future extension would be to have a model 
that takes into account the vectorial nature of the wave field. 

The power spectral density of the wave field is illustrated in Figure ya and can be qualified 
as being continuous and essentially featureless. Notice that all frequencies are expressed in the 



spacecraft reference frame, in which they are Doppler-shifted by the strong solar wind. The spectral 
densities are almost the same for the two spacecraft. Figure ^3 shows the wave field probability 
distribution, which has non-Gaussian tails. The departure from Gaussianity should be underlined, 
since it is a necessary condition for having nonlinear wave-wave interactions [Kim and Powers, 
1979]. 




normalized amplitude 

Figure 3: (a) Power spectral density and (b) probability distribution of the magnetic field of 
AMPTE-UKS projected along the direction of maximum variance. The solid line corresponds to a 
Gaussian distribution with the same variance; all amplitudes are normalized to have unit variance. 



3 Modeling the Nonlinear Transfer Function 

Much work has been done on the theory of nonlinear transfer functions in turbulence [e.g., Monin 
and Yaglom, 1975; Krommes, 1997] but relatively little is known about their inference from exper- 
imental data, which can be an unwieldy task. Early results were obtained in the context of neutral 
fluid turbulence [Uberoi, 1963; Van Atta and Chen, 1969; Lii et al., 1982; Ritz et al., 1988a] and 
later in plasmas [Ritz and Powers, 1986; Ritz et al., 1988b; Ritz et al., 1989; Kim et al., 1996]. 
Powers, Ritz and their coworkers contributed to the development of a computational framework for 
two-point measurements [Ritz and Powers, 1986; Ritz et al., 1989], thereby rendering the technique 



easily accessible to a large class of experiments. Their results, however, have remained overlooked, 
presumably because of the apparent computational investment and the difficulty in validating es- 
timates that are prone to errors. In this paper we show how to increase the robustness of the 
estimates by using continuous wavelet transforms instead of the usual Fourier transform. 

3.1 Volterra Series 

Consider a stationary wave field which is measured in time and in space, and let y{x,t) denote 
fluctuations around a fixed value. We are interested in describing the dynamics of this wave field 
with the following general model: 

^ = n.(-.')). (1) 

where F(y) is a continuous, nonlinear and time-invariant operator. Wiener [1958] showed that for 
a large class of causal systems F(y) can be expanded as a Volterra (or Volterra- Frechet- Wiener) 
series [Schetzen, 1980], which we write here after taking the Fourier transform of the time variable 



dyix,uj) 
dx 



^T{uj)y{x,u) (2) 

T{uji,uj2) y{x,uji)y{x,uj2) 
X S(uJl + LU2 — Lo) dujiduj2 
T{uJi,uj2, LOs) y{x, ui) y{x, UJ2) 
X y{x, 073) S(uji + LU2 + <-^3 — Lo) duJiduj2dL03, 



The kernels F are directly related to the higher-order spectra of the process and have a physical 
meaning. r(w), T{uji,uj2)j and T{u!i, 102,103) are respectively called the linear, quadratic, and cubic 
interaction terms. The generic situation corresponds to a leading linear term, which describes the 
linear dynamics of the system, such as the linear growth rate and the dispersion. The quadratic 
term expresses three- wave processes in which interactions occur within triads of waves that satisfy 
the resonance condition 

CO = oji + L02 ■ (3) 

The cubic term similarly describes four-wave processes whose frequencies satisfy the selection rules 

LO = LOl + LO2 + ^3 or LO + LOl = LO2 + ^3 ■ (4) 

The main motivation for using a Volterra series expansion stems from its ability to describe var- 
ious weakly nonlinear processes in plasmas [Kadomtsev, 1982], ranging from generic drift wave 
turbulence [Balk et al., 1990; Horton and Hasegawa, 1994] to Langmuir turbulence as described by 
the Zakharov equations [Musher et al., 1995]. Particular attention has been given to Hamiltonian 
systems [Zakharov et al, 1985] in which the kernels can be calculated explicitly. The resonant 



interactions defined by (||) and (^ are further known to be the building elements of turbulence as 
observed in collisionless plasmas: the decay and modulational instabilities, for example, are ade- 
quately described in terms of three- wave and four- wave interactions [Krasnosel'skikh and Lefeuvre, 
1993]. 

Theoretical and experimental considerations show that for weak turbulence the low-order 
Volterra kernels are the predominant ones. Indeed, the characteristic timescale associated with 
the action of a gth-order kernel increases with q, making low-order kernels much more likely to 
rule the dynamics [Zakharov et al., 1985]. In practice, (g) may thus safely be truncated after the 
cubic term and quite often even a quadratically nonlinear model suffices. 

3.2 Strong Versus Weak Turbulence 

The nonlinear model of (|2|) formally applies to weak turbulence only, in which the dispersion and 
the characteristic growth rates of the Fourier modes are small. Solar wind turbulence, on the other 
hand, is often considered as being of the strong turbulence type. The region we study is actually a 
mixture between the two since the dynamical properties of the wave field are dominated by a small 
population of energetic ions interacting with a plasma of the weak turbulence type. The weakness 
of the dispersion [Dwrfofc de Wit et al.,, 1995] and the relatively small value of the linear growth 
rate (see Section g|) support the validity of the weak turbulence approximation here. 

The extension from weak to strong turbulence as a first approximation implies a loosening of 
the resonance conditions (equations (H)-(||)) to account for the finite bandwidth of the wave packets 
[Horton and Hasegawa, 1994]. We shall take this spectral broadening implicitly into account by 
projecting the wave field on wavelets instead of Fourier modes (see Section ^. 

3.3 Spatial Versus Temporal Description 

Equation (yj) is actually a particular case of a class of models that describe both the spatial 
and the temporal structure of the wave field. In a more general setting, both wavenumbers and 
frequencies must satisfy resonance conditions. Energy and momentum conservation force three- 
wave interactions to occur along the resonant manifold 

w(ki-Fk2)=CJl(ki)+CJ2(k2) , (5) 

where k denotes the wave-number vector. 

The wavenumber dependence of the interaction is often omitted by lack of spatially resolved 
experiments. It raises an important point, which is the separation between spatial and temporal 
scales, and the distinction between stationarity and homogeneity. In our experiment, the wave field 
is convected past the satellites by the solar wind and so the angular frequency ujsat we observe in 
the spacecraft frame is in fact Doppler-shifted, giving, tOsat = i-^pi + k • "Vgw However, because the 
wavenumber k and the solar wind velocity Vsw are almost parallel, and because of the fast solar 
wind, we may approximate tOsat ~ '^pi + k Vsw ~ k Vgw ■ This expression suggests that the spatial 
structure of the wave field, projected on the solar wind velocity vector, can be probed simply by 
measuring the time evolution and viceversa. This approximation is known as the Taylor hypothesis 



and allows us to exchange temporal dynamics and spatial structure 



dt 



dy_ 
' dx 



(6) 



There now remains to convert the Eulerian representation of the experiment into a Lagrangian 
one, in which the magnetic structures are followed from one spacecraft to the other. A space-time 
representation of the spacecraft (see Figure Q) indeed reveals that by taking the difference of the 
spacecraft signals, we mix the wave field time derivative and the spatial gradient. A Galilean 
transformation is needed y{x, t) -^ y{x, t' = t — x/vsw), which we do by shifting the AMPTE-IRM 
time scries by r = —5x/vsw — —0.39 s. An additional correction of —0.28 s is needed to compensate 
for differences in timing conventions [see Schwartz et al., 1992]. 
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Figure 4: Representation of the spacecraft in space-time, showing the correspondence between 
spatial separation Sx and time delay r. The actual position of IRM must be moved to IRM* to 
compensate for the effect of the solar wind. 



3.4 Inferring the Model Coefficients 

Physical insight into our model can be gained by introducing the real-valued density of waves 
E(u;,x), by analogy to the number of quasi particles in condensed matter theory. For a large 
ensemble of waves with different frequencies, the random-phase approximation holds, giving 



{y{uji,x)y*{uj2,x)) = E{ll>i,x)S{uji - UJ2) 



(7) 



Angle brackets denote ensemble-averaging, which is often replaced by time averaging, assuming 
ergodicity. From equations (y) and ([7|) we obtain the kinetic equation 



dEiuj, x) „ , s, ^, X 



dt 



(8) 



T{uJi,U!2) S{u!i + UJ2 — Uj) duji duj2 



which models the nonhnear evolution taking place between the two spacecraft. Notice that we 
used the Taylor hypothesis (equation (||)) to interchange spatial and temporal derivatives. The 
quantities of interest are the average linear growth rate in time 

7(c^) = v,^ Re [r{Lu)] (9) 

and the average quadratic energy transfer rate 

T{uJi,U;2) ^VsniRe[r{uJi,UJ2) (10) 

X {y{uji,x)y{uj2,x)y*{uji+uj2,x))] . 

The latter attests the existence of nonlocal interactions which are a hallmark of nonlinearity. 
Equation (||) shows that the energy flux in Fourier space dE/dt results from a balance between 
energy dissipation (or gain) at a given frequency lu and spectral energy transfers between lu and 
other frequencies. 

Nonlinear transfer functions have the advantage of revealing both the magnitude and the orien- 
tation of spectral energy fluxes: positive values of T(a;i, W2) correspond to three-wave interactions 
in which spectral components with angular frequencies wi > and 0^2 > transfer energy to 
the component uj — uji + uj2- We shall write this a.s uJi + uj2 ^ ^- Conversely, negative values 
correspond to decay processes uj ^ uii + uj2. 

There is a close resemblance between the definition of the energy transfer function (equation 
(p!o|)) and that of the auto bispectrum [Mendel, 1991], 

B{uji,uj2) ^ {y{uji,x) y{uj2,x) y*{uji + uj2,x)) , (11) 

which has been widely used for quantifying quadratic wave interactions in plasmas [Kim and 
Powers, 1979; Lagoutte et al., 1989; LaBelle and Lund, 1992; Pecseli et at, 1993; Dudok de Wit 
and Krasnosel'skikh, 1995; Bale et al., 1996]. Transfer functions, however, are more informative 
since they detect the presence of nonlinear interactions between the observation points irrespective 
of what happened farther upstream. Consider for example a wave field that underwent nonlinear 
interactions during its early history but is now fully static. This wave field will have a nonzero cross- 
bispectrum even though nonlinear interactions aren't actually taking place. A nonzero bispectrum 
thus does not necessarily attest the existence of wave-wave interactions at the observation point. 
Such a caveat was put forward in an analytic example by Pecseli and Trulsen [1993]. The energy 
transfer function on the contrary detects whether energy is being exchanged between spectral 
modes, causing the amplitudes and the phases to vary locally in time and in space. This new 
information can be accessed only by comparing the wave field as it goes from one observation 
point to the other. Finally, we note that the existence of energy transfers presupposes a weak 
nonstationarity or inhomogeneity of the wave field. 

3.5 Symmetries 

The real-valued nature of the data and the definition of the transfer function automatically give 
rise to a number of symmetry relations, which shrink the principal domain in frequency space, 
significantly reducing the number of model coefficients to be computed [see Nam and Powers, 



1994]. The principal domain of the energy transfer function is even smaller, since for uji + L02 = uj 
we have 

T{loi,uJ2) = T{lu2, LOi) = -T{uj, -LOl) = -T {uj , -CJ2) • (12) 

4 Estimating the Nonlinear Transfer Function 

Our principal problem is the robust estimation of Volterra kernels from finite and noise-corrupted 
data. This problem may be alleviated by assuming a Gaussian probability distribution of the wave 
field, since the different regressors can then be identified separately. This assumption, however, 
rarely holds in practice. Incidentally, it is precisely the nonlinearity that causes the distribution to 
depart from Gaussianity. We therefore follow a more general procedure along the line developed 
by Ritz and Powers [1986] and later improved by Kim and Powers, [1988]. 

For discrete values of the frequency and with two-point measurements, (0) becomes 

yu,{x+5x) - y^{x) 

7 = i^u,yuj{x) (13) 

+ 2 ^ ^^i,'^2 Vu^ii^) Vi^Ax) 

ij = ui+u2 

+ ... 

where yuj{x) is the discrete Fourier transform of y{x,t) and {m} = {aJi,(xi2, • • • ji^jv^} are the 
regularly spaced frequencies. Without loss of generality we assume that the ensemble average 
vanishes {y{x,t)) = 0. It is convenient to express the complex wave field yu{x) as 

y^{x) = |y^(x)|e^"^-(^) . (14) 

From (|l3| ) and ( |l4|) and in the limit where (5a; ^ 0, we obtain a new system [Ritz and Powers, 
1986] 

Y^ = L^U^ + \ Y, Q",,^2t^-it^-2 + --- (15) 



with 





uj^uj^ +^2 


Y^ -~ 


= yuj{x + Sx) 


u^ -- 


= yUx) 


L^ -- 


= {T^St + 1 ^ J Sq^^) , 


WuJi,U2 


= K u Ste^^'f'-^ 

UJl .UJ2 


6t = 


= Sx/Vsw 


S(f>ui = 


= 't'ujix + Sx) - (t>^{x) 



(16) 



From this system, the physical quantities F^ and A'^_^ ^^ can be computed directly, as shown by 
Ritz et al. [1989]. 



Equation (051) formally represents a nonlinear transfer function that links an output Y^ (the 
waveform of AMPTE-IRM) to an input U^ (the waveform of AMPTE-UKS). The estimation of 
the linear part of such a transfer function is a central problem in system identification, for which 
well-established techniques exist [Ljung, 1987; Priestley, 1981]. Comparatively few experimental 
efforts, however, have been directed toward the robust estimation of quadratic and higher-order 
transfer functions [Tick, 1961; Brillinger, 1970; Billings, 1980; Bendat, 1990]. For the sake of 
simplicity, we shall henceforth restrict ourselves to quadratically nonlinear models. 

The simplest solution consists in selecting the model whose coefficients minimize the squared 
residual errors e^i between the measured wave field and the predicted one Y^ 



j^ iji j^ 111 



(17) 



The problem then reduces to a multiple linear regression with a unique solution. For each angular 
frequency uj we solve for H^j, 

(18) 



UojHu — Yi^ 



with 



U. 






t/c.i(2)C/„-„i(2) 






(19) 



Numbers refer to different ensembles collected under identical conditions, T denotes transposition 
and the number of unknown coefficients is Nc- The conceptual simplicity of this approach and its 
straightforward generalization to cubic and higher-order interactions are clear advantages. 

The nonlinear transfer function can be obtained by solving the overdetermined set of equations 
(equation (|l^)) using conventional least squares techniques but deeper insight can be gained by 
multiplying these equations on the left by U* , giving 



{\U. 



{u^ u:, u: 



{u: u^' u^-. 



([/;, u:_^, u^., u^.^.,) 



{U*Y^ 









(20) 



The leading matrix (also called higher-order autocovariance matrix) can be divided into four 
blocks, one with a second order moment (the power spectral density), two with third-order moments 
(the bispectra) and one with fourth-order moments. The fact that moments of various orders are 
needed to properly estimate the linear properties of the wave field recalls the well-known closure 
problem which is ubiquitous in the spectral modelling of turbulence. Equation (20) also shows 
how the non-Gaussian nature of the wave field enters the results. If the wave field were Gaussian, 
then the off-diagonal blocks of the higher-order autocovariance matrix would vanish and a separate 
estimation of the different Volterra kernels would be possible. 
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5 Wavelet Versus Fourier Transform 

For the solution of (20) to be numerically stable and physically relevant, it is essential to have 
^ens 3> Nc {Nc is the number of unknown coefhcients) and U^^ nonsingular. A compromise is 
thus needed between N^ns and the number of different Fourier modes N^, hereafter referred to 
as the mesh resolution. Usually, time series are divided into (possibly overlapping) sequences, 
each of which is Fourier transformed. A better compromise can be achieved with wavelets, which 
offer additional resolution in time at the expense of a lower frequency resolution. The continuous 
wavelet transform of y{t) is defined as 



I'^'^Ta'*{^)'' 



y{a,T)^ j y{t)-^h*^-—jdt, (21) 

where h{t) is the analyzing wavelet and a its scale. The optimum tradeoff between time and 
frequency resolution is achieved with Gaussian or Morlet wavelets 

hit) ^ ^^TTiW^'"' ^"^'"^ ' ^'') 

for which each scale is related to an instantaneous angular frequency uj — 271 /a. The frequency 
resolution, defined in terms of the cutoff frequency at 3 dB is Alo/uj = 1/4ct and the usual Fourier 
transform is recovered for cr ^ oo. 

Compared to windowed Fourier transforms, the wavelet transforms yield statistically better 
behaved estimates of the spectral properties [van Milligen et al., 1995; Dudok de Wit and Kras- 
nosel'skikh, 1995]. Our motivation, however, is not just computational but also stems from the 
ability of wavelets to resolve transient and soliton-like features [Farge et al., 1996]. Indeed, the 
strongly turbulent magnetic field shows transient structures that are more akin to wavelets than 
to coherent waves with an infinite extension. 

The main drawback of this approach is its greater computational burden, since the number of 
ensembles N^^s now almost equals the number of samples. Furthermore, we are left with a free 
parameter, the wavelet width a. Since a fixed mesh resolution is wanted, with no spectral overlap 
between adjacent components 1^. and VL +i, we adapt the wavelet width to the frequency in order 
to have a > lo /ASuo. An additional condition cr > 1 is imposed to prevent the analyzing wavelet 
from being too much distorted. 

6 Validation Criteria for the Transfer Function Model 

Validation is a key issue in Volterra model identification. There exists no single satisfactory cri- 
terion for performing such a validation, but to a large extent we can rely on well-established 
techniques that have been developed for linear systems [Ljung, 1987]. 

Since our problem involves the solution of a linear system of equations, a good starting point 
is an inspection of the degree of independence between the columns of the matrix U^j (equation 
(19)) and the output Y^^. The correlation function between Y^^ and the first column of U;^ 

11 



indicates how well the linear transfer function succeeds in predicting the output. This is the 
coherence function, which is bounded between and 1. Likewise, the correlation function between 
the output and other columns of the matrix gives 

This is the cross-bicoherence, i.e., the cross- bispectrum normalized to the power spectral density. 
Its value is bounded between zero for uncorrelated waves and unity for triads of waves whose 
phases are totally correlated. For a cubic model, one similarly defines the cross-tricoherence, 
which quantifies the strength of four-wave interactions 

7c('^1.^2,W3) = -— |2W|r. rj rj |2N " (^S) 

Note that there exist variants of these definitions [e.g. Kravtchenko-Berejnoi et al., 1995]. 

Higher-order coherence functions reveal which spectral components arc likely to be involved 
in nonlinear interactions. They do not, however, tells us whether the model is actually good in 
predicting the wave field. A high bicoherence, for example, does not yet justify the choice of 
a quadratically nonlinear model. A more global figure of merit is obtained by comparing the 
measured output signal Y^ to the model prediction Y^ 



2 ,,,^_ \{Y^y:)\ 



iivi^) - ,^::i ,.. • (26) 



In practice, one half of the data is used to estimate the transfer function while the other half is 
kept for cross-validation. These simple prescription tools can be complemented by tests on the 
residuals, etc. 

7 Choosing the Right Model 

The choice of the model parameters involves three main issues. More specific aspects are deferred 
to the appendix. 

7.1 Choice of the Model Order 

The basic question of the model order should ideally be answered by computing Volterra kernels 
for various orders and truncating the series as soon as they become negligible. The finite size of 
the data does not allow this, so the question should rather be: How faithfully does a truncated 
low-order model reproduce the observed dynamics ? 

As mentioned before, there are several reasons to believe that a low-order model should capture 
most of the dynamics, especially in weak turbulence. This can be verified in different ways. Ritz 
and Powers [1986] considered nonlinear correlations between linear and nonlinear terms. We focus 
instead on the predictive capacity of the models, using the cross-validation defined in d2q). We 
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built first-, second-, and third-order models, all of which were tested against the data (the third 
order model could could not have as much frequency resolution because of its large number of 
degrees of freedom) . 

Figure shows the result of the cross-validation applied to a linear and to a quadratically 
nonlinear model. Both models succeed relatively well in predicting the low frequency part of 
the AMPTE-IRM waveform. The performance drop with increasing frequency is a well-known 
effect, which cannot be compensated simply by increasing the model complexity. Possible causes 
are the decreasing signal-to-noise ratio, the finite lifetime and the dispersion of the wave packets, 
fluctuations in the solar wind velocity, and the 1-D approximation of our model. 

The central result here is the close performance of the linear and the quadratic models, which 
attests the predominantly linear behavior of the wave field and a priori supports the choice of a 
low-order model. A notable exception occurs around / = w/27r « 0.5 Hz, where a quadratic model 
brings some improvement. This shall see later that nonlinear effects are indeed important in that 
frequency band. 



7.2 Choice of the Frequency Range 

There is a strong impetus for reducing as much as possible the number of degrees of freedom 
of our model. One way of doing this is by reducing the frequency range. As shown in Figure ^ 
fluctuations with frequencies beyond 0.8 Hz cannot be satisfactorily modeled and so one may safely 
truncate the frequency range at 1 Hz, above which the power spectral content becomes negligible 
anyway. We checked that higher-order coherence functions vanish as well above 1 Hz. 

A further reduction in the number of degrees of freedom can in principle be achieved by dis- 
carding in the linear system (equation (|l8[)) those columns of the matrix which are not significantly 
correlated with the output Y^ . Such a reduction is permitted when the nonlinear interactions are 
very localized in frequency space. 



nonlinear model 
linear model 




0.5 1 1.5 2 

frequency [Hz] 

Figure 5: Cross-coherence between the measured and the simulated output, for a linear and a 
nonlinear model. Values close to or below the bias level are not considered to be significant 
[Bendat and Piersol, 1986]. 



13 



7.3 Choice of the Mesh Resolution 

The frequency spacing 5u) (which is proportional to l/N^,) must be small enough to distinguish 
important features such as spectral lines and yet as large as possible to prevent the model from 
being overdetermined. Since we deal with broadband turbulence, a relatively coarse mesh should 
a priori suffice. Nonlinear parametric models [Billings, 1980] may be needed when closely spaced 
lines must be resolved. 




20 40 60 

mesh resolution 

Figure 6: Condition number of the matrix Uuj versus (a) frequency with a fixed mesh resolution 
Ni^ = 20, and (b) versus N^ for a fixed frequency / « 0.5 Hz. 



The impact of the mesh resolution is best revealed by the condition nmnbcr [Golub and Van 
Loan, 1993] of the matrix U^^, which gives a figure of merit for the ill-posedness of ([l8|). The 
condition number is at best 1 and typically should not exceed a few hundreds; its value is displayed 
in Figure g for different mesh resolutions. The condition degrades for increasing N^^ because 
more coefficients have to be estimated from the same sample; another reason is the increasing 
coUinearity between columns of Ui^. These problems may be partially alleviated by projecting U(^ 
on an orthogonal basis [see Im et al., 1996]. 

From these considerations, we choose a quadratically nonlinear model with N^ = 40 and a 
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frequency range from to 1 Hz. 

8 Linear Properties of the Wave Field 

We now focus on the interpretation of the model and start with the leading term, which is the 
linear one. The linear kernel T,^ can conveniently be split into a real and an imaginary part 



r{L0)^-fiu)+j9{uj) 



(27) 



The imaginary part 



0H«-^im[iog(t/:y^)] 



(28) 



expresses the average phase-shift undergone by the wave-packets as they move from one spacecraft 
to the other. It is therefore related to the wavenumber k, averaged over the power spectral density, 
hyk^9/vsw. 
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Figure 7: (a) Imaginary and (b) real parts of the linear transfer function, as calculated using a 
linear and a nonlinear model. Error bars correspond to ± one standard deviation. The results 
are inaccurate above 0.6 Hz because the model fails to reproduce small-amplitude fluctuations 
correctly. 
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The average dispersion relation k{uj) is shown in Figure u[ Given our timing convention, we 
are in the plasma rest frame and positive wavenumbers correspond to a sunward motion. Error 
bars correspond to one standard deviation as calculated from the least-squares fit of the model. 

Figure ^ shows that the wave field is essentially dispersionless up to about 0.5 Hz. Above this 
frequency, the dispersion becomes positive, and high-frequency waves move ahead of low-frequency 
ones. We refer to previous work [Dudok de Wit et al., 1995] for a discussion on this, but just note 
that the modeling fails above 0.6 Hz, presumably because of the low power spectral density. 

The real part of the Volterra kernel 

gives the linear growth averaged as before over k. The results are shown in Figure M. The negative 
value of 7(w) attests a damping of the waves, so we conclude that the wave field is on average 
linearly stable. An exception occurs below 0.2 Hz, where the wave field grows as it goes from one 
spacecraft to the other. Although this growth rate is subject to a rather large uncertainty, its 
positive sign is statistically significant. Interestingly, this unstable frequency band coincides with 
that of the SLAMS and therefore lends strong support to the instability of these structures. The 
unstable nature of the SLAMS has already been conjectured [Schwartz and Burgess, 1991], but we 
now have the first direct evidence for a dynamic evolution. 

From the linear growth rate, one can estimate the characteristic time needed for the SLAMS to 
grow, assuming that there is no nonlinear mechanism to saturate such a growth. We find t = I/7 w 
10 s, a value that should be compared to the characteristic lifetime of these structures 

This ratio is sufficiently small to justify a linearization of the growth process (and the weak turbu- 
lence approximation) and yet large enough to make the instability of the SLAMS easily detectable. 
It is instructive to check here the assumption of nonlinearity by comparing these results to 
what we would obtain by fitting a linear model. The two growth rates are compared in Figure M 
and a discrepancy appears. We attribute this to the energy redistribution process between Fourier 
modes, which is neglected in the linear model and correctly taken into account in the nonlinear 
one. As we shall see later, the linear model tries to compensate an energy flux around 0.5 Hz by 
artificially lowering the damping rate. The use of a nonlinear model for assessing linear properties 
should therefore not be underestimated. 

9 Nonlinear Properties of the Wave Field 

We now consider the properties of the second-order Volterra kernel. As mentioned before, this is 
the only nonlinear term we can reliably estimate given the amount of data. 
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9.1 Phase Couplings 



The second-order Volterra kernels r(a;i,a;2) as such are not very informative, even though they 
contain all the relevant information on the quadratic couplings. Better insight can be gained by 
looking at the the cross-bicoherence (equation (|24|)), which is indicative of the strength of the 
quadratic interactions. In the same way, we compute the cross- tricoherence (equation (|2^)) to 
study cubic interactions, even though the third-order kernel itself cannot be reliably estimated. 
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Figure 8: The cross-bicoherence, displayed in the principal domain, for frequency-adding interac- 
tions only (/i, /2 > 0). Its value is bounded between and 1. 
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Figure 9: Cross-tricoherence, shown here for interactions of the type /i + /2 + /3 ~ / Hz only, with 
/ii/2,/3 > 0. Values are bounded between and 1, and the dashed line refers to the principal 
domain. 
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The cross-bicoherence and cross-tricoherence are displayed in Figures @ and respectively. In 
both Figures ^ and || the support is restricted to the nonredundant and positive frequency domain. 
The most conspicuous result is the presence of local maxima that attest the existence of phase 
couplings between specific spectral modes. We conclude from the cross-bicoherence that a signif- 
icant phase coupling occurs between wave packets whose frequencies satisfy the summation rule 
0.1-1- 0.45 = 0.55 Hz. The cross-tricoherence reveals a significant coupling for 0.1 + fi + f^ = 0.55 
Hz, with 0.1 < fi < fm < 0.55 Hz. Both couphngs relate wave packets whose frequencies are about 
0.1 Hz and 0.55 Hz, with possibly some intermediate frequencies to enable the phase coupling. As 
shown in previous work [Dudok de Wit and KrasnoseVskikh, 1995], these characteristic frequencies 
respectively correspond to the SLAMS and to the discrete whistler wave packets that frequently 
occur at the leading edge of SLAMS. The nonzero cross-bicoherence and cross-tricoherence thus 
reveal the existence of a causal relationship between the SLAMS and the whistlers. 

9.2 Energy Transfers 

The last step now consists in determining whether the SLAMS and the whistlers are exchanging 
energy or if they are just remnants of a process that took place farther upstream. To do so, we 
compute the quadratic energy transfer function, shown in Figure HO. A significant energy flux 
appears at 0.1 -f 0.45 -^ 0.55 Hz, which corresponds to an energy transfer going from the SLAMS 
to the whistlers. This is the central result of our paper, from which we conclude that the whistlers 
are much more likely to be a decay product of the SLAMS than some instability triggered by 
them. Such a conclusion was partly anticipated, but only the energy transfer function can give 
unambiguous evidence for it. 
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Figure 10: Spectral energy transfer rate (in arbitrary units) using the same representation as for 
the cross-bicoherence. The confidence interval is about equal to the spacing between two contour 
levels. 



The other patterns in Figure 10 also have an interpretation. The 0.1 + 0.1 -^ 0.2 Hz transfer 



corresponds a first harmonic generation due the nonlinear steepening of the SLAMS. 
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9.3 Power Balance 

Further insight into the dynamics of the wave field can be gained by studying the power balance. 
Consider the truncated second-order wave kinetic equation (equation m)) 



dt 



^l.E^ 



E ^-^ 



(31) 






Stationarity [dE^/dt = 0) is approximately reached when the two terms on the right-hand side 
cancel; these two terms are plotted in Figure |n|. 
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Figure 11: Relative change of spectral power showing the contributions of the linear growth term 
(27^), and of the three-wave interactions {2'Y^T^^_^^/E^). The results are unreliable above 0.6 
Hz. 

At low frequencies (/ < 0.2 Hz) the growth rate and the energy transfer indeed approximately 
cancel each other and so the wave field amplitude should not vary much in time. We conclude 
that the decay of the SLAMS is approximately compensated by their linear instability. At higher 
frequencies the power balance becomes increasingly negative, suggesting that the waves are on 
average damped. A plausible damping mechanism would be resonant particle damping, but the 
various approximations made in our model may actually cause the damping to be overestimated 
at high frequencies (see the appendix). 

9.4 Interpretation 

A coherent scenario now emerges, which is schematized in Figure n3. The SLAMS appear as 
dynamically evolving structures that progressively grow out of the wave field by drawing energy 
from energetic ions. As they grow, nonlinear effects enter into play. The transfer function analysis 
shows that the dominant process is a nonlinear wave interaction that compensates the growth by an 
energy transfer toward high-frequency whistler waves (some energy may also go into low-frequency 
waves). The whistler waves in turn move ahead of the SLAMS because of the positive dispersion 
and are eventually damped by dissipation. 
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Figure 12: Schematic representation of the power spectral density of the magnetic field, showing 
where the energy enters the wave field and where it is transferred before being dissipated. The 
SLAMS are located in the hatched zone. 

The emergence of such right-handed circularly polarized waves out of the left-handed linearly 
polarized SLAMS shows some striking similarities with the expected behavior of solitary waves 
[Hada et at, 1989] and also recalls the behavior of shock fronts in weakly dispersive media [Karp- 
man, 1975]. All these phenomena have in common a competition between dispersion and nonlin- 
earity, whose distinctive manifestation is the resilience of the shape of the SLAMS. 

To finish, let us visualize how the nonlinear interactions show up in the time domain. Figure |l3| 
represents a particular SLAMS with the measured and the predicted wave field. We decomposed 
the latter into its linear and quadratic contributions. As expected, the linear contribution captures 
most of the dynamics but does not correctly reproduce the fast oscillations at the leading edge of 
the SLAMS and which correspond to a distorted whistler wave packet. A quadratic contribution 
is definitely needed here to fit the observations. 

10 Conclusions 

This study reveals how Volterra models can be used to infer nonlinear properties from a turbulent 
wave field using two-point measurements. Provided the model is carefully validated, it can give 
direct access to the wave field growth rate and to the energy transfer function. 

We used the Volterra approach to analyze plasma turbulence as observed just upstream the 
Earth's quasi-parallel bow shock by the AMPTE spacecraft. An important feature of the wave 
field is the occurrence of nonlinear magnetic structures termed SLAMS. Our analysis attests the 
coexistence of two competing mechanisms: the SLAMS progressively grow by drawing energy 
from hot ions but before overturning they saturate and release the excess of energy into high 
frequency whistler waves that move ahead of them due to dispersion. The dynamical evolution of 
the SLAMS and their differential velocity [Schwartz et al., 1992] support the conjecture in which 
they supposedly merge into an extended front that constitutes the bow shock. 

The method we advocate here is applicable to other types of events provided they are recorded 
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Figure 13: Excerpt of Figure 2, showing one component of the magnetic field with a typical SLAMS 
at t =70-76 s and its whistler precursor at t =76-82 s. (a) The wave field measured by AMPTE- 
IRM is compared with the model prediction, (b) The model prediction is split into its linear and 
quadratic constituents. 

by multiple-spacecraft missions whose configuration satisfies some constraints (see the appendix). 
A number of improvements can be made, such as the generalization to vector fields, which would 
allow anisotropy effects to be included. In some cases the addition of a source term that enforces 
the stationarity of the wave field may be desirable. 
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A Appendix: Limitations of the Technique 

The nonlinear transfer function cannot be meaningfully assessed without keeping in mind several 
limitations and potential pitfalls. The most important ones are listed here. 
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A.l One-Dimensional Approach 

Restricting the analysis to two probes means that we can only study structures propagating par- 
allel to the spacecraft separation vector. Structures propagating obliquely to it will cause an 
overestimation of the damping rate and an underestimation of the energy transfer. 

A. 2 Anisotropy 

Our scalar model can in principle be generalized to vectorial quantities with the interesting per- 
spective of addressing anisotropy effects. The price to pay for this is a larger number of degrees 
of freedom. Neglecting the vectorial nature of the variables may actually alter the power balance 
because of the omission of the wave field rotation. Although the rotation we measure is quite small, 
we believe this effect to be sufficient to accentuate the globally negative trend observed in Figure 

A. 3 Taylor Hypothesis 

Another problem comes from the connection between frequency and wavenumber space, for which 
the Taylor hypothesis was invoked at the beginning. The dispersionless approximation remains 
valid up to about 0.5 Hz. Above this limit, the dispersion becomes positive and hence the resonance 
conditions (equations (H)-(0)) are altered. The finite frequency resolution of the wavelets can easily 
accommodate such small detunings in the resonance conditions. 

A. 4 Vahdity of the Model 

The weakest point of our approach is the difficulty in justifying the validity of a low-order model 
in a definite way. A quadratically nonlinear model suffices for describing the main features of 
the wave field (section |7|), but Figure y reminds us that cubic interactions cannot be neglected. 
Although we are confident in the conclusions drawn from our second-order model, one should keep 
in mind that the results remain approximate as long as cubic and possibly even higher-order effects 
are not taken into account. 

A. 5 Calibration of the Probes 

It is essential that the two probes (the magnetometers here) be properly calibrated and have the 
same instrumental transfer function for our analysis to be meaningful. Although this is not a 
problem for the frequency range we are considering, it may exclude diagnostics that have either 
an unreliable calibration or a nonlinear response. 
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